tic
clear;
load 'C:\Users\scott\Dropbox\Spatial Probit PSRM\PSRM_R&R\Monte_Code\Results_Final\sp_tscs_phi50_rho25';

X = store_x0(nt-(n-1):nt,1);

clearvars -except X; 

load 'C:\Users\scott\Dropbox\Spatial Probit PSRM\PSRM_R&R\Monte_Code\Results_Final\sp_tscs_phi30_rho25';
clear std; 

X_prod = ones(t,1); 
X_steady = X*X_prod';
X_steady = reshape(X_steady,nt,1);


time = 10; 

%Immediate Effects 
unit_1 = 51; 
unit_2 = 72; 

%pt = pt';        

trials = 100; 

mean_diff_within = zeros(time,trials);
sdev_diff_within = zeros(time,trials);

for h=1:trials

        
n_freq = 100;
draws_fx = mvnrnd(pt0(:,h)', vcov_hat0(:,:,h), n_freq);

cp_diff = zeros(time,n_freq); 

for k=1:n_freq

mult_fx = inv(i_nt - draws_fx(k,3)*WNT - draws_fx(k,4)*TL);
SIGMA = mult_fx*mult_fx';
rf_cut_full = mult_fx*(draws_fx(k,1) + draws_fx(k,2)*X_steady); 

p_1 = zeros(time,1);
p0_1 = normcdf(rf_cut_full(unit_2),0,sqrt(SIGMA(unit_2,unit_2)));

for i=1:time
        p_1(i) = mvncdf([-Inf,-Inf], [rf_cut_full([unit_2]), rf_cut_full([n*(i-1) + unit_1])],0,SIGMA([unit_2,n*(i-1) + unit_1],[unit_2,n*(i-1) + unit_1])); 
end

cp_1 = zeros(time,1);


for i = 1:time
    cp_1(i) = p_1(i)/p0_1;
end

p_0 = zeros(time,1);
p0_0 = normcdf(-rf_cut_full(unit_2),0,sqrt(SIGMA(unit_2,unit_2)));

for i=1:time
        p_0(i) = mvncdf([rf_cut_full([unit_2]),-Inf], [Inf, rf_cut_full([n*(i-1) + unit_1])],0,SIGMA([unit_2,n*(i-1) + unit_1],[unit_2,n*(i-1) + unit_1])); 
end

cp_0 = zeros(time,1);


for i = 1:time
    cp_0(i) = p_0(i)/p0_0;
end
cp_diff(:,k) = cp_1-cp_0;

end

mean_diff_within(:,h) = mean(cp_diff'); 
sdev_diff_within(:,h) = std(cp_diff');
end

mean_diff_across = mean(mean_diff_within'); 
sdev_diff_across = std(mean_diff_within'); 
se_diff_across = mean(sdev_diff_within'); 

toc

save 'C:\Users\scott\Dropbox\Spatial Probit PSRM\PSRM_R&R\Monte_Code\Results_Final\sp_tscs_phi30_rho25_fx_results'; 